# import counts
countsData <- read.delim(file = "../01_input/all.counts", sep = " ")
# preview counts
head(countsData)
##                  chr start stop strand length embryo_cells_rep1
## WBGene00014450 MtDNA     1   55      +     55                 0
## WBGene00014451 MtDNA    58  111      +     54                 0
## WBGene00010957 MtDNA   113  549      +    437                 0
## WBGene00010958 MtDNA   549  783      +    235                 0
## WBGene00014452 MtDNA   785  840      +     56                 0
## WBGene00014453 MtDNA   842  896      +     55                 0
##                embryo_cells_rep2 embryo_GFPminus_rep1 embryo_GFPminus_rep2
## WBGene00014450                 0                    0                    0
## WBGene00014451                 0                    0                    0
## WBGene00010957                 0                    0                    0
## WBGene00010958                 0                    0                    0
## WBGene00014452                 0                    0                    0
## WBGene00014453                 0                    0                    0
##                embryo_GFPminus_rep3 embryo_GFPplus_rep1 embryo_GFPplus_rep2
## WBGene00014450                    0                   0                   0
## WBGene00014451                    0                   0                   0
## WBGene00010957                    0                   0                   0
## WBGene00010958                    0                   0                   0
## WBGene00014452                    0                   0                   0
## WBGene00014453                    0                   0                   0
##                embryo_GFPplus_rep3 embryo_whole_rep2 embryo_whole_rep3
## WBGene00014450                   0                 0                 0
## WBGene00014451                   0                 0                 0
## WBGene00010957                   0                 0                 0
## WBGene00010958                   0                 0                 0
## WBGene00014452                   0                 0                 0
## WBGene00014453                   0                 0                 0
##                L1_cells_rep1 L1_cells_rep2 L1_cells_rep3 L1_GFPminus_rep1
## WBGene00014450             0             0             0                0
## WBGene00014451             0             0             0                0
## WBGene00010957             0             0             0                0
## WBGene00010958             0             0             0                0
## WBGene00014452             0             0             0                0
## WBGene00014453             0             0             0                0
##                L1_GFPminus_rep2 L1_GFPminus_rep3 L1_GFPplus_rep1
## WBGene00014450                0                0               0
## WBGene00014451                0                0               0
## WBGene00010957                0                0               0
## WBGene00010958                0                0               0
## WBGene00014452                0                0               0
## WBGene00014453                0                0               0
##                L1_GFPplus_rep2 L1_GFPplus_rep3 L1_whole_rep1 L1_whole_rep2
## WBGene00014450               0               0             0             0
## WBGene00014451               0               0             0             0
## WBGene00010957               0               0             0             0
## WBGene00010958               0               0             0             0
## WBGene00014452               0               0             0             0
## WBGene00014453               0               0             0             0
##                L1_whole_rep3 L3_cells_rep1 L3_cells_rep2 L3_cells_rep3
## WBGene00014450             0             0             0             0
## WBGene00014451             0             0             0             0
## WBGene00010957             0             0             0             0
## WBGene00010958             0             0             0             0
## WBGene00014452             0             0             0             0
## WBGene00014453             0             0             0             0
##                L3_GFPminus_rep1 L3_GFPplus_rep2 L3_GFPminus_rep3
## WBGene00014450                0               0                0
## WBGene00014451                0               0                0
## WBGene00010957                0               0                0
## WBGene00010958                0               0                0
## WBGene00014452                0               0                0
## WBGene00014453                0               0                0
##                L3_GFPplus_rep1 L3_GFPminus_rep2 L3_GFPplus_rep3 L3_whole_rep1
## WBGene00014450               0                0               0             0
## WBGene00014451               0                0               0             0
## WBGene00010957               0                0               0             0
## WBGene00010958               0                0               0             0
## WBGene00014452               0                0               0             0
## WBGene00014453               0                0               0             0
##                L3_whole_rep2 L3_whole_rep3
## WBGene00014450             0             0
## WBGene00014451             0             0
## WBGene00010957             0             0
## WBGene00010958             0             0
## WBGene00014452             0             0
## WBGene00014453             0             0
# print samples
colnames(countsData[6:ncol(countsData)])
##  [1] "embryo_cells_rep1"    "embryo_cells_rep2"    "embryo_GFPminus_rep1"
##  [4] "embryo_GFPminus_rep2" "embryo_GFPminus_rep3" "embryo_GFPplus_rep1" 
##  [7] "embryo_GFPplus_rep2"  "embryo_GFPplus_rep3"  "embryo_whole_rep2"   
## [10] "embryo_whole_rep3"    "L1_cells_rep1"        "L1_cells_rep2"       
## [13] "L1_cells_rep3"        "L1_GFPminus_rep1"     "L1_GFPminus_rep2"    
## [16] "L1_GFPminus_rep3"     "L1_GFPplus_rep1"      "L1_GFPplus_rep2"     
## [19] "L1_GFPplus_rep3"      "L1_whole_rep1"        "L1_whole_rep2"       
## [22] "L1_whole_rep3"        "L3_cells_rep1"        "L3_cells_rep2"       
## [25] "L3_cells_rep3"        "L3_GFPminus_rep1"     "L3_GFPplus_rep2"     
## [28] "L3_GFPminus_rep3"     "L3_GFPplus_rep1"      "L3_GFPminus_rep2"    
## [31] "L3_GFPplus_rep3"      "L3_whole_rep1"        "L3_whole_rep2"       
## [34] "L3_whole_rep3"
# import metadata and process file
metadata1 <- read.table(file = "../01_input/RWP27_metadata.tsv", header = FALSE, stringsAsFactors = FALSE) %>% bind_rows(read.table(file = "../01_input/RWP26_metadata.tsv", header = FALSE, stringsAsFactors = FALSE)) %>%
  bind_rows(read.table(file = "../01_input/RWP30_metadata.tsv", header = FALSE, stringsAsFactors = FALSE))

colnames(metadata1) <- c("Filename.Fwd", "Filename.Rev", "names")
head(metadata1)
##           Filename.Fwd         Filename.Rev                names
## 1 RW57_S10_L003_R1_001 RW57_S10_L003_R2_001    embryo_cells_rep1
## 2 RW58_S11_L003_R1_001 RW58_S11_L003_R2_001  embryo_GFPplus_rep1
## 3 RW59_S12_L003_R1_001 RW59_S12_L003_R2_001 embryo_GFPminus_rep1
## 4 RW60_S13_L003_R1_001 RW60_S13_L003_R2_001    embryo_whole_rep2
## 5 RW61_S14_L003_R1_001 RW61_S14_L003_R2_001    embryo_cells_rep2
## 6 RW62_S15_L003_R1_001 RW62_S15_L003_R2_001  embryo_GFPplus_rep2
# separate and process sample info
metadata1 <- metadata1 %>% separate(names, sep = "_", into = c("stage", "sample", "rep"), remove = FALSE)
metadata1 <- metadata1 %>% mutate(stage = fct_relevel(stage, c("embryo", "L1", "L3")), 
                     sample = fct_relevel(sample, c("whole", "cells", "GFPplus", "GFPminus")),
                     rep = fct_relevel(rep, c("rep1", "rep2", "rep3")),
                     names = fct_relevel(names, metadata1$names)
                     )

# Order columns according to metadata1 order
countsData <- countsData  %>% select(chr:length, sort(metadata1$names))
head(countsData)
##                  chr start stop strand length embryo_cells_rep1
## WBGene00014450 MtDNA     1   55      +     55                 0
## WBGene00014451 MtDNA    58  111      +     54                 0
## WBGene00010957 MtDNA   113  549      +    437                 0
## WBGene00010958 MtDNA   549  783      +    235                 0
## WBGene00014452 MtDNA   785  840      +     56                 0
## WBGene00014453 MtDNA   842  896      +     55                 0
##                embryo_GFPplus_rep1 embryo_GFPminus_rep1 embryo_whole_rep2
## WBGene00014450                   0                    0                 0
## WBGene00014451                   0                    0                 0
## WBGene00010957                   0                    0                 0
## WBGene00010958                   0                    0                 0
## WBGene00014452                   0                    0                 0
## WBGene00014453                   0                    0                 0
##                embryo_cells_rep2 embryo_GFPplus_rep2 embryo_GFPminus_rep2
## WBGene00014450                 0                   0                    0
## WBGene00014451                 0                   0                    0
## WBGene00010957                 0                   0                    0
## WBGene00010958                 0                   0                    0
## WBGene00014452                 0                   0                    0
## WBGene00014453                 0                   0                    0
##                embryo_whole_rep3 embryo_GFPplus_rep3 embryo_GFPminus_rep3
## WBGene00014450                 0                   0                    0
## WBGene00014451                 0                   0                    0
## WBGene00010957                 0                   0                    0
## WBGene00010958                 0                   0                    0
## WBGene00014452                 0                   0                    0
## WBGene00014453                 0                   0                    0
##                L1_whole_rep1 L1_cells_rep1 L1_GFPplus_rep1 L1_GFPminus_rep1
## WBGene00014450             0             0               0                0
## WBGene00014451             0             0               0                0
## WBGene00010957             0             0               0                0
## WBGene00010958             0             0               0                0
## WBGene00014452             0             0               0                0
## WBGene00014453             0             0               0                0
##                L1_whole_rep2 L1_cells_rep2 L1_GFPplus_rep2 L1_GFPminus_rep2
## WBGene00014450             0             0               0                0
## WBGene00014451             0             0               0                0
## WBGene00010957             0             0               0                0
## WBGene00010958             0             0               0                0
## WBGene00014452             0             0               0                0
## WBGene00014453             0             0               0                0
##                L1_whole_rep3 L1_cells_rep3 L1_GFPplus_rep3 L1_GFPminus_rep3
## WBGene00014450             0             0               0                0
## WBGene00014451             0             0               0                0
## WBGene00010957             0             0               0                0
## WBGene00010958             0             0               0                0
## WBGene00014452             0             0               0                0
## WBGene00014453             0             0               0                0
##                L3_whole_rep1 L3_cells_rep1 L3_GFPplus_rep1 L3_GFPminus_rep1
## WBGene00014450             0             0               0                0
## WBGene00014451             0             0               0                0
## WBGene00010957             0             0               0                0
## WBGene00010958             0             0               0                0
## WBGene00014452             0             0               0                0
## WBGene00014453             0             0               0                0
##                L3_whole_rep2 L3_cells_rep2 L3_GFPminus_rep2 L3_GFPplus_rep2
## WBGene00014450             0             0                0               0
## WBGene00014451             0             0                0               0
## WBGene00010957             0             0                0               0
## WBGene00010958             0             0                0               0
## WBGene00014452             0             0                0               0
## WBGene00014453             0             0                0               0
##                L3_whole_rep3 L3_cells_rep3 L3_GFPplus_rep3 L3_GFPminus_rep3
## WBGene00014450             0             0               0                0
## WBGene00014451             0             0               0                0
## WBGene00010957             0             0               0                0
## WBGene00010958             0             0               0                0
## WBGene00014452             0             0               0                0
## WBGene00014453             0             0               0                0
# Generate a table called "cts" out of the countsData table. Subset the countsData.
cts <- as.matrix(countsData %>% select(metadata1$names))
head(cts)
##                embryo_cells_rep1 embryo_GFPplus_rep1 embryo_GFPminus_rep1
## WBGene00014450                 0                   0                    0
## WBGene00014451                 0                   0                    0
## WBGene00010957                 0                   0                    0
## WBGene00010958                 0                   0                    0
## WBGene00014452                 0                   0                    0
## WBGene00014453                 0                   0                    0
##                embryo_whole_rep2 embryo_cells_rep2 embryo_GFPplus_rep2
## WBGene00014450                 0                 0                   0
## WBGene00014451                 0                 0                   0
## WBGene00010957                 0                 0                   0
## WBGene00010958                 0                 0                   0
## WBGene00014452                 0                 0                   0
## WBGene00014453                 0                 0                   0
##                embryo_GFPminus_rep2 embryo_whole_rep3 embryo_GFPplus_rep3
## WBGene00014450                    0                 0                   0
## WBGene00014451                    0                 0                   0
## WBGene00010957                    0                 0                   0
## WBGene00010958                    0                 0                   0
## WBGene00014452                    0                 0                   0
## WBGene00014453                    0                 0                   0
##                embryo_GFPminus_rep3 L1_whole_rep1 L1_cells_rep1 L1_GFPplus_rep1
## WBGene00014450                    0             0             0               0
## WBGene00014451                    0             0             0               0
## WBGene00010957                    0             0             0               0
## WBGene00010958                    0             0             0               0
## WBGene00014452                    0             0             0               0
## WBGene00014453                    0             0             0               0
##                L1_GFPminus_rep1 L1_whole_rep2 L1_cells_rep2 L1_GFPplus_rep2
## WBGene00014450                0             0             0               0
## WBGene00014451                0             0             0               0
## WBGene00010957                0             0             0               0
## WBGene00010958                0             0             0               0
## WBGene00014452                0             0             0               0
## WBGene00014453                0             0             0               0
##                L1_GFPminus_rep2 L1_whole_rep3 L1_cells_rep3 L1_GFPplus_rep3
## WBGene00014450                0             0             0               0
## WBGene00014451                0             0             0               0
## WBGene00010957                0             0             0               0
## WBGene00010958                0             0             0               0
## WBGene00014452                0             0             0               0
## WBGene00014453                0             0             0               0
##                L1_GFPminus_rep3 L3_whole_rep1 L3_cells_rep1 L3_GFPplus_rep1
## WBGene00014450                0             0             0               0
## WBGene00014451                0             0             0               0
## WBGene00010957                0             0             0               0
## WBGene00010958                0             0             0               0
## WBGene00014452                0             0             0               0
## WBGene00014453                0             0             0               0
##                L3_GFPminus_rep1 L3_whole_rep2 L3_cells_rep2 L3_GFPminus_rep2
## WBGene00014450                0             0             0                0
## WBGene00014451                0             0             0                0
## WBGene00010957                0             0             0                0
## WBGene00010958                0             0             0                0
## WBGene00014452                0             0             0                0
## WBGene00014453                0             0             0                0
##                L3_GFPplus_rep2 L3_whole_rep3 L3_cells_rep3 L3_GFPplus_rep3
## WBGene00014450               0             0             0               0
## WBGene00014451               0             0             0               0
## WBGene00010957               0             0             0               0
## WBGene00010958               0             0             0               0
## WBGene00014452               0             0             0               0
## WBGene00014453               0             0             0               0
##                L3_GFPminus_rep3
## WBGene00014450                0
## WBGene00014451                0
## WBGene00010957                0
## WBGene00010958                0
## WBGene00014452                0
## WBGene00014453                0
# Reorganize the metadata table so the names2 column are now headers
rownames(metadata1)<- metadata1$names
coldata <- metadata1[,c("names", "stage", "sample", "rep")]
rownames(coldata) <- as.vector(metadata1$names)
# make grouping variable
coldata$group <- factor(paste0(coldata$stage, coldata$sample))
coldata
##                                     names  stage   sample  rep          group
## embryo_cells_rep1       embryo_cells_rep1 embryo    cells rep1    embryocells
## embryo_GFPplus_rep1   embryo_GFPplus_rep1 embryo  GFPplus rep1  embryoGFPplus
## embryo_GFPminus_rep1 embryo_GFPminus_rep1 embryo GFPminus rep1 embryoGFPminus
## embryo_whole_rep2       embryo_whole_rep2 embryo    whole rep2    embryowhole
## embryo_cells_rep2       embryo_cells_rep2 embryo    cells rep2    embryocells
## embryo_GFPplus_rep2   embryo_GFPplus_rep2 embryo  GFPplus rep2  embryoGFPplus
## embryo_GFPminus_rep2 embryo_GFPminus_rep2 embryo GFPminus rep2 embryoGFPminus
## embryo_whole_rep3       embryo_whole_rep3 embryo    whole rep3    embryowhole
## embryo_GFPplus_rep3   embryo_GFPplus_rep3 embryo  GFPplus rep3  embryoGFPplus
## embryo_GFPminus_rep3 embryo_GFPminus_rep3 embryo GFPminus rep3 embryoGFPminus
## L1_whole_rep1               L1_whole_rep1     L1    whole rep1        L1whole
## L1_cells_rep1               L1_cells_rep1     L1    cells rep1        L1cells
## L1_GFPplus_rep1           L1_GFPplus_rep1     L1  GFPplus rep1      L1GFPplus
## L1_GFPminus_rep1         L1_GFPminus_rep1     L1 GFPminus rep1     L1GFPminus
## L1_whole_rep2               L1_whole_rep2     L1    whole rep2        L1whole
## L1_cells_rep2               L1_cells_rep2     L1    cells rep2        L1cells
## L1_GFPplus_rep2           L1_GFPplus_rep2     L1  GFPplus rep2      L1GFPplus
## L1_GFPminus_rep2         L1_GFPminus_rep2     L1 GFPminus rep2     L1GFPminus
## L1_whole_rep3               L1_whole_rep3     L1    whole rep3        L1whole
## L1_cells_rep3               L1_cells_rep3     L1    cells rep3        L1cells
## L1_GFPplus_rep3           L1_GFPplus_rep3     L1  GFPplus rep3      L1GFPplus
## L1_GFPminus_rep3         L1_GFPminus_rep3     L1 GFPminus rep3     L1GFPminus
## L3_whole_rep1               L3_whole_rep1     L3    whole rep1        L3whole
## L3_cells_rep1               L3_cells_rep1     L3    cells rep1        L3cells
## L3_GFPplus_rep1           L3_GFPplus_rep1     L3  GFPplus rep1      L3GFPplus
## L3_GFPminus_rep1         L3_GFPminus_rep1     L3 GFPminus rep1     L3GFPminus
## L3_whole_rep2               L3_whole_rep2     L3    whole rep2        L3whole
## L3_cells_rep2               L3_cells_rep2     L3    cells rep2        L3cells
## L3_GFPminus_rep2         L3_GFPminus_rep2     L3 GFPminus rep2     L3GFPminus
## L3_GFPplus_rep2           L3_GFPplus_rep2     L3  GFPplus rep2      L3GFPplus
## L3_whole_rep3               L3_whole_rep3     L3    whole rep3        L3whole
## L3_cells_rep3               L3_cells_rep3     L3    cells rep3        L3cells
## L3_GFPplus_rep3           L3_GFPplus_rep3     L3  GFPplus rep3      L3GFPplus
## L3_GFPminus_rep3         L3_GFPminus_rep3     L3 GFPminus rep3     L3GFPminus
# Check that the names match  --> Should be TRUE
all(rownames(coldata) == colnames(cts))
## [1] TRUE

Make DESeqDataSet

Generate the DESeqDataSet. The variables in this design formula will be the type of sample, and the preparation date. This should reduce the variability between the samples based on when they were made.

From the vignette: “In order to benefit from the default settings of the package, you should put the variable of interest at the end of the formula and make sure the control level is the first level.”

The variable of interest is the sample type.

Using DESeqDataSetFromMatrix since I used the program featureCounts.

dds <- DESeqDataSetFromMatrix(countData = cts,
                              colData = coldata,
                              design = ~ group)

Visualize read count distribution

hist(log(rowSums(counts(dds))))
abline(v = log(10), col = "red", lty = 2)

Filter genes with low read counts

keep <- rowSums(counts(dds)) >=10
dds <- dds[keep,]
dds
## class: DESeqDataSet 
## dim: 26557 34 
## metadata(1): version
## assays(1): counts
## rownames(26557): WBGene00021406 WBGene00021407 ... WBGene00199694
##   WBGene00044951
## rowData names(0):
## colnames(34): embryo_cells_rep1 embryo_GFPplus_rep1 ... L3_GFPplus_rep3
##   L3_GFPminus_rep3
## colData names(5): names stage sample rep group

Perform Differential Expression

dds <- DESeq(dds)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
resultsNames(dds)
##  [1] "Intercept"                           "group_embryoGFPminus_vs_embryocells"
##  [3] "group_embryoGFPplus_vs_embryocells"  "group_embryowhole_vs_embryocells"   
##  [5] "group_L1cells_vs_embryocells"        "group_L1GFPminus_vs_embryocells"    
##  [7] "group_L1GFPplus_vs_embryocells"      "group_L1whole_vs_embryocells"       
##  [9] "group_L3cells_vs_embryocells"        "group_L3GFPminus_vs_embryocells"    
## [11] "group_L3GFPplus_vs_embryocells"      "group_L3whole_vs_embryocells"

Sample-to-sample distance matrix

vsd <- vst(dds, blind = FALSE)
sampleDists <- dist(t(assay(vsd)))
sampleDistMatrix <- as.matrix(sampleDists)
rownames(sampleDistMatrix) <- vsd$names
colnames(sampleDistMatrix) <- NULL
colors <- colorRampPalette( rev(brewer.pal(9, "Blues")) )(255)
pheatmap(sampleDistMatrix,
         clustering_distance_rows = sampleDists,
         clustering_distance_cols = sampleDists,
         col = colors)

vsd.corr.per.stage <- function(x, main){
vsd <- assay(vsd)[,metadata1 %>% filter(grepl(x, names)) %>% pull(names)]
sampleDists <- dist(t(vsd))
sampleDistMatrix <- as.matrix(sampleDists)
colnames(sampleDistMatrix) <- NULL
colors <- colorRampPalette( rev(brewer.pal(9, "Blues")) )(255)
pheatmap(sampleDistMatrix,
         clustering_distance_rows = sampleDists,
         clustering_distance_cols = sampleDists,
         col = colors, 
         main = main)
}
vsd.corr.per.stage("embryo", "Embryo Stage Intestine FACS RNA-seq Correlation")

vsd.corr.per.stage("L1", "L1 Stage Intestine FACS RNA-seq Correlation")

vsd.corr.per.stage("L3", "L3 Stage Intestine FACS RNA-seq Correlation")

All embryo stage pairwise comparisons

Code testing, proliminary visualization

#for loop attempt one
combinations<- data.frame(combos = 1:4, top = c("embryowhole", "embryoGFPplus", "embryoGFPminus", "embryoGFPplus"), bottom = c("embryocells", "embryocells", "embryocells", "embryoGFPminus"))

res_embryowhole_vs_embryocells <- results(dds, contrast = c("group", combinations$top[1], combinations$bottom[1]))
res_embryowhole_vs_embryocells_shrunk <- lfcShrink(dds, contrast = c("group", combinations$top[1], combinations$bottom[1]), type = "ashr", quiet = TRUE)

res_embryowhole_vs_embryocells <- results(dds, contrast = c("group", combinations$top[1], combinations$bottom[1]))
res_embryowhole_vs_embryocells_shrunk <- lfcShrink(dds, contrast = c("group", combinations$top[1], combinations$bottom[1]), type = "ashr", quiet = TRUE)

for(i in 1:nrow(combinations)){
  res.var.name<- paste("res",combinations$top[i], "vs", combinations$bottom[i], sep = "_")
  res.shrunk.var.name <- paste("res",combinations$top[i], "vs", combinations$bottom[i], "shrunk", sep = "_")
  print(res.var.name)
  print(res.shrunk.var.name)
  assign(res.var.name, results(dds, contrast = c("group", combinations$top[i], combinations$bottom[i])))
}
## [1] "res_embryowhole_vs_embryocells"
## [1] "res_embryowhole_vs_embryocells_shrunk"
## [1] "res_embryoGFPplus_vs_embryocells"
## [1] "res_embryoGFPplus_vs_embryocells_shrunk"
## [1] "res_embryoGFPminus_vs_embryocells"
## [1] "res_embryoGFPminus_vs_embryocells_shrunk"
## [1] "res_embryoGFPplus_vs_embryoGFPminus"
## [1] "res_embryoGFPplus_vs_embryoGFPminus_shrunk"
# res_embryoGFPplus_vs_embryoGFPminus
# plotMA(res_embryowhole_vs_embryocells_shrunk, ylim = ylim, main = "whole vs cells")
# combinations$top[1]
# class(res_embryowhole_vs_embryocells)

# for loop attempt two
samples <- paste("embryo", c("whole", "GFPplus", "cells", "GFPminus"), sep = "")
combos <- combn(samples, 2, simplify = FALSE)

for(i in 1:length(combos)){
  print(c(combos[[i]][1], combos[[i]][2]))
  res.var.name<- paste("res",combos[[i]][1], "vs", combos[[i]][2], sep = "_")
  res.shrunk.var.name <- paste("res",combos[[i]][1], "vs", combos[[i]][2], "shrunk", sep = "_")
  print(res.var.name)
  print(res.shrunk.var.name)
  assign(res.var.name, results(dds, contrast = c("group", combos[[i]][1], combos[[i]][2])))
  assign(res.shrunk.var.name, lfcShrink(dds, contrast = c("group", combos[[i]][1], combos[[i]][2]), type = "ashr", quiet = TRUE))
}
## [1] "embryowhole"   "embryoGFPplus"
## [1] "res_embryowhole_vs_embryoGFPplus"
## [1] "res_embryowhole_vs_embryoGFPplus_shrunk"
## [1] "embryowhole" "embryocells"
## [1] "res_embryowhole_vs_embryocells"
## [1] "res_embryowhole_vs_embryocells_shrunk"
## [1] "embryowhole"    "embryoGFPminus"
## [1] "res_embryowhole_vs_embryoGFPminus"
## [1] "res_embryowhole_vs_embryoGFPminus_shrunk"
## [1] "embryoGFPplus" "embryocells"  
## [1] "res_embryoGFPplus_vs_embryocells"
## [1] "res_embryoGFPplus_vs_embryocells_shrunk"
## [1] "embryoGFPplus"  "embryoGFPminus"
## [1] "res_embryoGFPplus_vs_embryoGFPminus"
## [1] "res_embryoGFPplus_vs_embryoGFPminus_shrunk"
## [1] "embryocells"    "embryoGFPminus"
## [1] "res_embryocells_vs_embryoGFPminus"
## [1] "res_embryocells_vs_embryoGFPminus_shrunk"
combos
## [[1]]
## [1] "embryowhole"   "embryoGFPplus"
## 
## [[2]]
## [1] "embryowhole" "embryocells"
## 
## [[3]]
## [1] "embryowhole"    "embryoGFPminus"
## 
## [[4]]
## [1] "embryoGFPplus" "embryocells"  
## 
## [[5]]
## [1] "embryoGFPplus"  "embryoGFPminus"
## 
## [[6]]
## [1] "embryocells"    "embryoGFPminus"
par(mfrow = c(3, 3), mar=c(2,2,1,1))
ylim = c(-10,10)
plotMA(res_embryowhole_vs_embryocells_shrunk, ylim = ylim, main = "whole vs cells")
plotMA(res_embryowhole_vs_embryoGFPplus_shrunk, ylim = ylim, main = "whole vs GFP+")
plotMA(res_embryowhole_vs_embryoGFPminus_shrunk, ylim = ylim, main = "whole vs GFP-")
plotMA(res_embryoGFPplus_vs_embryocells_shrunk, ylim = ylim, main = "GFP+ vs cells")
plotMA(res_embryoGFPplus_vs_embryoGFPminus_shrunk, ylim = ylim, main = "GFP+ vs GFP-")
plotMA(res_embryocells_vs_embryoGFPminus_shrunk, ylim = ylim, main = "cells vs GFP-")

ggplot version of plotting

all_embryo_comparisons <- data.frame()
for(i in 1:length(combos)){
tobind <- as.data.frame(lfcShrink(dds, contrast = c("group", combos[[i]][1], combos[[i]][2]), type = "ashr", quiet = TRUE)) %>% 
  rownames_to_column(var = "WBGeneID") %>% 
  mutate(comparison = paste(combos[[i]][1], combos[[i]][2], sep = "_vs_")) 
all_embryo_comparisons <- bind_rows(all_embryo_comparisons, tobind)
}


ggplot(all_embryo_comparisons %>% filter(!is.na(padj)), aes(x = log10(baseMean), y = log2FoldChange, color = padj < 0.1)) +
  geom_point(shape = 16, alpha = 0.1, stroke = 0, size = 1) +
  ylim(c(-10,10))+
  facet_wrap(~comparison) +
  scale_color_manual(values = c("black", "red"), name = "q.value < 0.1") +
  theme_classic()
## Warning: Removed 19 rows containing missing values (geom_point).

# Functions

pairwise_array_df <- function(stage) {
  samples <- paste(stage, c("whole", "GFPplus", "cells", "GFPminus"), sep = "")
  combos <- combn(samples, 2, simplify = FALSE)
  # print(combos)
  all_pairwise_comparisons <- data.frame()
  for (i in 1:length(combos)) {
    tobind <-
      as.data.frame(lfcShrink(
        dds,
        contrast = c("group", combos[[i]][1], combos[[i]][2]),
        type = "ashr",
        quiet = TRUE
      )) %>%
      rownames_to_column(var = "WBGeneID") %>%
      mutate(comparison = paste(combos[[i]][1], combos[[i]][2], sep = "_vs_"))
    all_pairwise_comparisons <- bind_rows(all_pairwise_comparisons, tobind)
  }
  all_pairwise_comparisons
}

MA_plot_array <- function(in.df){
  ggplot(in.df %>% mutate(padj = replace_na(padj, 1)), aes(x = log10(baseMean), y = log2FoldChange, color = padj < 0.1)) +
  geom_point(shape = 16, alpha = 0.1, stroke = 0, size = 1) +
  ylim(c(-10,10))+
  facet_wrap(~comparison) +
  scale_color_manual(values = c("black", "red"), name = "q.value < 0.1") +
  theme_classic()
}

alt_hyp_res_df <- function(stage){
samples <- paste(stage, c("whole", "GFPplus", "cells", "GFPminus"), sep = "")
combos <- combn(samples, 2, simplify = FALSE)
thresh = 0.5
sig = 0.01
hyps = c("greater", "less", "lessAbs")
df <- data.frame()
for(i in 1:length(combos)){
  for(hyp in hyps){
    thresh_res <- results(dds, contrast = c("group", combos[[i]][1],combos[[i]][2]), lfcThreshold=thresh, altHypothesis = hyp, alpha = sig)
    tobind<-data.frame(plotMA(thresh_res, returnData = TRUE), 
                      comparison = paste(combos[[i]][1], combos[[i]][2], sep = "_vs_"), 
                       type = hyp)
    df <- bind_rows(df, tobind)
  }
}
df
}

de_category_MA_plot <- function(df){
df %>% filter(isDE == TRUE) %>%
  ggplot(aes(x = log10(mean), y = lfc, color = type)) +
  geom_point(data =df, shape = 16, alpha = 0.01, stroke = 0, size = 1, color = "grey") +
  geom_point(shape = 16, alpha = 0.5, stroke = 0, size = 1) +
  ylim(c(-10,10))+
  facet_wrap(~comparison) +
  # scale_color_manual(values = c("black", "red"), name = "q.value < 0.1") +
  theme_classic()
}

options(dplyr.summarise.inform = FALSE)
de_category_bar_plot <- function(df){
  df %>% filter(isDE == TRUE) %>% group_by(comparison, type) %>% summarize(genes = n()) %>%
  ggplot(aes(x = type, y = genes, label = genes, fill = type)) +
  geom_bar(stat = "identity") +
  geom_text(vjust = -0.25) +
  facet_wrap(~comparison) +
    theme_classic()
}
embryo_alt_hyp_res_df <- alt_hyp_res_df("embryo")
de_category_MA_plot(embryo_alt_hyp_res_df)
## Warning: Removed 195 rows containing missing values (geom_point).
## Warning: Removed 44 rows containing missing values (geom_point).

de_category_bar_plot(embryo_alt_hyp_res_df)

L1_alt_hyp_res_df<- alt_hyp_res_df("L1")
de_category_MA_plot(L1_alt_hyp_res_df)
## Warning: Removed 12 rows containing missing values (geom_point).

de_category_bar_plot(L1_alt_hyp_res_df)

L3_alt_hyp_res_df<- alt_hyp_res_df("L3")
de_category_MA_plot(L3_alt_hyp_res_df)
## Warning: Removed 135 rows containing missing values (geom_point).
## Warning: Removed 27 rows containing missing values (geom_point).

de_category_bar_plot(L3_alt_hyp_res_df)

Embryo Intestine FACS MA plots

embryo_pairwise_res_shrunk <- pairwise_array_df(stage = "embryo")

MA_plot_array(embryo_pairwise_res_shrunk)
## Warning: Removed 21 rows containing missing values (geom_point).

L1 Intestine FACS MA plots

L1_pairwise_res_shrunk <- pairwise_array_df(stage = "L1")

MA_plot_array(L1_pairwise_res_shrunk)
## Warning: Removed 3 rows containing missing values (geom_point).

L3 Intestine FACS MA plots

L3_pairwise_res_shrunk <- pairwise_array_df(stage = "L3")
MA_plot_array(L3_pairwise_res_shrunk)
## Warning: Removed 33 rows containing missing values (geom_point).

DESeq2 builtin plotting function

res_embryoGFPplus_vs_embryoGFPminus <- results(dds, contrast = c("group", "embryoGFPplus", "embryoGFPminus"))
res_L1GFPplus_vs_L1_GFPminus <- results(dds, contrast = c("group", "L1GFPplus", "L1GFPminus"))
res_L3GFPplus_vs_L3_GFPminus <- results(dds, contrast = c("group", "L3GFPplus", "L3GFPminus"))
par(mfrow=c(1,3),mar=c(2,2,1,1))
ylim <- c(-15,15)
# drawLines <- function() abline(h=c(-2,2),col="dodgerblue",lwd=2)
sig = 0.01
plotMA(res_embryoGFPplus_vs_embryoGFPminus, ylim=ylim, main = "Embryo GFP+ vs GFP-", alpha = sig)
plotMA(res_L1GFPplus_vs_L1_GFPminus, ylim=ylim, main = "L1 GFP+ vs GFP-", alpha = sig)
plotMA(res_L3GFPplus_vs_L3_GFPminus, ylim=ylim, main = "L3 GFP+ vs GFP-", alpha = sig)

res_embryoGFPplus_vs_embryoGFPminus_apeglm <- lfcShrink(dds, contrast = c("group", "embryoGFPplus", "embryoGFPminus"), type = "ashr")
## using 'ashr' for LFC shrinkage. If used in published research, please cite:
##     Stephens, M. (2016) False discovery rates: a new deal. Biostatistics, 18:2.
##     https://doi.org/10.1093/biostatistics/kxw041
res_L1GFPplus_vs_L1GFPminus_apeglm <- lfcShrink(dds, contrast = c("group", "L1GFPplus", "L1GFPminus"), type = "ashr")
## using 'ashr' for LFC shrinkage. If used in published research, please cite:
##     Stephens, M. (2016) False discovery rates: a new deal. Biostatistics, 18:2.
##     https://doi.org/10.1093/biostatistics/kxw041
res_L3GFPplus_vs_L3GFPminus_apeglm <- lfcShrink(dds, contrast = c("group", "L3GFPplus", "L3GFPminus"), type = "ashr")
## using 'ashr' for LFC shrinkage. If used in published research, please cite:
##     Stephens, M. (2016) False discovery rates: a new deal. Biostatistics, 18:2.
##     https://doi.org/10.1093/biostatistics/kxw041
par(mfrow=c(1,3),mar=c(2,2,1,1))
ylim <- c(-10,10)
# drawLines <- function() abline(h=c(-2,2),col="dodgerblue",lwd=2)
sig = 0.01
plotMA(res_embryoGFPplus_vs_embryoGFPminus_apeglm, ylim=ylim, main = "Embryo GFP+ vs GFP-", alpha = sig)
plotMA(res_L1GFPplus_vs_L1GFPminus_apeglm, ylim=ylim, main = "L1 GFP+ vs GFP-", alpha = sig)
plotMA(res_L3GFPplus_vs_L3GFPminus_apeglm, ylim=ylim, main = "L3 GFP+ vs GFP-", alpha = sig)

plotCounts(dds, "WBGene00001578", intgroup = "group", main = "ges-1 read counts")

plotCounts(dds, "WBGene00001578", intgroup = "group", returnData = TRUE) %>% 
  separate(group, sep = "(?<=embryo)|(?<=L1)|(?<=L3)", into = c("stage", "sample"), remove = FALSE) %>%
  ggplot(aes(x = sample, y = count)) +
  geom_boxplot() +
  geom_point() +
  facet_grid(~stage) +
  ggtitle("ges-1 read counts") +
  theme_classic()

plotCounts(dds, "WBGene00001250", intgroup = "group", main = "elt-2 read counts")

plotCounts(dds, "WBGene00001250", intgroup = "group", returnData = TRUE) %>% 
  separate(group, sep = "(?<=embryo)|(?<=L1)|(?<=L3)", into = c("stage", "sample"), remove = FALSE) %>%
  ggplot(aes(x = sample, y = count)) +
  geom_boxplot() +
  geom_point() +
  facet_grid(~stage) +
  ggtitle("elt-2 read counts") +
  theme_classic()

# Annotate and quantify tissue specific genes

tissue_specific_genes <- read_csv(file = "../../01_tissue_specific_genes/03_output/tissue_specific_genes_220202.csv", show_col_types = FALSE)


tissue_annotated_MA <- function(in_res){
df <- as.data.frame(in_res) %>% rownames_to_column(var = "WBGeneID") %>%
  left_join(tissue_specific_genes, by = "WBGeneID") %>%
  mutate(tissue = replace_na(tissue, "other"))

df %>%
  ggplot(aes(x = log10(baseMean), y = log2FoldChange, color = tissue)) +
  geom_point(data =df %>% select(-tissue), shape = 16, alpha = 0.1, stroke = 0, size = 1, color = "grey") +
  geom_point(shape = 16, alpha = 0.5, stroke = 0, size = 1) +
  facet_wrap(~tissue) +
  ylim(c(-10,10)) +
  theme_classic()
}
tissue_annotated_MA(res_embryoGFPplus_vs_embryoGFPminus_shrunk)
## Warning: Removed 96 rows containing missing values (geom_point).
## Warning: Removed 12 rows containing missing values (geom_point).

tissue_annotated_MA(res_L1GFPplus_vs_L1GFPminus_apeglm)
## Warning: Removed 8 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_point).

tissue_annotated_MA(res_L3GFPplus_vs_L3GFPminus_apeglm)
## Warning: Removed 48 rows containing missing values (geom_point).
## Warning: Removed 6 rows containing missing values (geom_point).

L3_apeglm_df <- as.data.frame(res_L3GFPplus_vs_L3GFPminus_apeglm) %>% rownames_to_column(var = "WBGeneID") %>%
  left_join(tissue_specific_genes, by = "WBGeneID") %>%
  mutate(tissue = replace_na(tissue, "other"), padj = replace_na(padj, 1))

L3_apeglm_df %>%
  ggplot(aes(x = log10(baseMean), y = log2FoldChange, color = (log2FoldChange>0.5 & padj < 0.01)|(log2FoldChange<0.5 & padj < 0.01))) +
  # geom_point(data =df %>% select(-tissue), shape = 16, alpha = 0.1, stroke = 0, size = 1, color = "grey") +
  geom_point(shape = 16, alpha = 0.5, stroke = 0, size = 1) +
  # facet_wrap(~tissue) +
  ylim(c(-10,10)) +
  theme_classic()
## Warning: Removed 6 rows containing missing values (geom_point).

tissue_gene_quant <- function(res, sig = 0.01, thresh = 0.5){
  my_plot <- as.data.frame(res) %>% rownames_to_column(var = "WBGeneID") %>%
  left_join(tissue_specific_genes, by = "WBGeneID") %>%
  mutate(tissue = replace_na(tissue, "other"), padj = replace_na(padj, 1)) %>% mutate(status = case_when(
  log2FoldChange > thresh & padj < sig ~ "greater",
  log2FoldChange < thresh & padj < sig ~ "less",
  TRUE ~ "no_sig_diff"
)) %>%
  group_by(tissue, status) %>%
  summarise(genes = n()) %>%
  ggplot(aes(x = status, y = genes, label = genes)) +
  geom_bar(stat = "identity") +
  geom_text(vjust = -0.25) +
  facet_wrap(~tissue)+
  ggtitle(deparse(substitute(res)))
  my_plot
}
tissue_gene_quant(res_embryoGFPplus_vs_embryoGFPminus)

tissue_gene_quant(res_L1GFPplus_vs_L1_GFPminus)

tissue_gene_quant(res_L3GFPplus_vs_L3_GFPminus)